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We present an efficient perturbative method to obtain both static and dynamic polarizabilities 
and hyperpolarizabilities of complex electronic systems. This approach is based on the solution 
of a frequency dependent Sternheimer equation, within the formalism of time-dependent density 
functional theory, and allows the calculation of the response both in resonance and out of resonance. 
Furthermore, the excellent scaling with the number of atoms opens the way to the investigation of 
response properties of very large molecular systems. To demonstrate the capabilities of this method, 
we implemented it in a real-space (basis-set free) code, and applied it to benchmark molecules, 
namely CO, H2O, and paranitroaniline (PNA). Our results are in agreement with experimental and 
previous theoretical studies, and fully validate our approach. 

PACS numbers: 



I. INTRODUCTION 

The optical properties of a material are essentially de- 
termined by the response of the electrons to an external 
field. If this applied field is small, the induced dipole of 
the system can be expanded in powers of the field.— ^ 
The first order coefficient is the so-called electric po- 
larizability a. This quantity describes, e.g., the dielec- 
tric properties of the material and how light is absorbed 
and emitted. In second order we obtain the first hy- 
perpolarizability (3, that is responsible for the processes 
of second-harmonic generation, optical rectification and 
Pockles effect* 1 - Higher order terms can be related to 
other electro-optical effects like third-harmonic genera- 
tion, the Kerr effect, etc. Finally, polarizabilities are 
called static or dynamic if the perturbing field is static 
or frequency dependent. 

The importance of the linear term, a, is well known in 
Physics and Chemistry i 3 - On the other hand, non-linear 
(i.e. beyond first order) optical effects have gained quite 
some interest lately due to their technological applica- 
tions in opto-electronic devices. Nonlinear optical ma- 
terials can be used to convert light to shorter (bluer) 
wavelengths, which can be focused to a smaller spot size. 
Shorter wavelength light sources would hence yield higher 
density optical recording media (such as DVDs and CDs) . 
Other applications include tunable light sources, image 
recognition systems and adaptive optics. 

Several methods have been used to calculate (hy- 
perpolarizabilities of finite systems In the static 
case, the simplest approach is finite differences^ that uses 
the definition of the polarizabilities as derivatives of the 
dipole (or of the total energy) with respect to the applied 
field. Calculations are performed at various (small) field 
strengths, and the required derivatives are evaluated nu- 
merically. This method is simple and straightforward to 



implement. However, it requires many total energy eval- 
uations, and these need to have a very high precision to 
obtain reasonable numerical derivatives. Moreover, it is 
not possible to generalize this idea to the dynamic case. 

Another widely used approach is perturbation theory, 
of which more than one flavor exists. In the sum-over- 
states method 4 -^ the (hyper)polarizabilities are written 
as an infinite sum over occupied and empty states, that 
involves the ground-state eigenvalues and dipole matrix 
elements. In a similar vein, one can obtain the polariza- 
bilities from the corresponding response functions writ- 
ten in the product basis of occupied and empty states 
or in terms of Green's functions^ Note that these tech- 
niques can be used both for static and dynamic response. 
Although widely used by the community, these meth- 
ods have several shortcomings. First, results are often 
difficult to converge with the number (and quality) of 
the empty states. Second, the scaling with the number 
of atoms is quite unfavorable, making hard the applica- 
tion to the study of large system, like nanostructures or 
molecules with biological interest. 

A different technique also used to obtain both static 
and dynamic linear polarizabilities is the direct solu- 
tion of the time-dependent Schrodinger equation in real 
time.— In this way only the occupied subspace is needed 
and the scaling with system size is excellent [0(N 2 ), 
where N is the number of atoms]. However, this method 
cannot be easily generalized to extract hyperpolarizabil- 
ities. 

A very interesting approach that essentially solves the 
problems mentioned above is the Sternheimer equation. 13 
Although a perturbative technique, it avoids the use of 
empty states, and has a quite good scaling with the 
number of atoms. This method has already been used 
for the calculation of many response properties^ like 
atomic vibrations (phonons), electron-phonon coupling, 
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magnetic response, etc. In the domain of optical re- 
sponse, this method has been mainly used for static re- 
sponse, although a few first principles calculations for low 
frequency (far from resonance) (hyper)polarizabilities 
have appeared] 14 i 15 ' 16 i 17 Recently, a reformulation of the 
Sternheimer equation in a super operator formalism was 
presented^ When combined with a Lanczos solver, it 
allows to calculate very efficiently the first order po- 
larizability for the whole frequency spectrum. However 
the generalization of this method to higher orders is not 
straight-forward. 

In this Article, we propose a modified version of the 
Sternheimer equation that is able to cope with both static 
and dynamic response in and out of resonance. The so- 
lution of the first-order Sternheimer equation gives us 
access to both a and (3. Higher order polarizabilities can 
be obtained from an hierarchy of Sternheimer equations. 
Exchange and correlation effects are treated at the level 
of density functional theory (DFT) 19 for static polari- 
zabilities and time-dependent DFT (TDDFT) 20 for the 
dynamic case. Compared to other quantum-chemistry 
approaches, density functional methods have a somewhat 
lower accuracy but are lighter numerically, allowing the 
study of much larger systems. In the present work we 
focus on finite systems but the method has also been 
applied to periodic systems for the non-resonant casei^I 
Note that in the field of DFT, the Sternheimer equa- 
tion is often referred to as density functional perturbation 
theory X 

The rest of this Article is organized as follows. In 
Sect. |H] we present the derivation of the frequency- 
dependent Sternheimer equation and show how to ob- 
tain the linear polarizability and first hyperpolarizability 
from its solution. In the following section we give some 
details concerning the implementation of our method. In 
Sect. |IV] we apply this theory to several test molecules, 
comparing our results to other calculations and experi- 
ments. Finally we present our conclusions and a brief 
outlook. 



II. THEORY 

A. Linear Response 

Within TDDFT, the quantum state of an interact- 
ing electronic system is described by the time-dependent 
Kohn-Sham equations (atomic units will be used unless 
explicitly stated) 

d 

i-^ifj m (r 7 t) = H KS (t)ip m (r 7 t) . (1) 

The Kohn-Sham Hamiltonian is written as 
V 2 

H K S = Y + W ext(f,*) +W H artree (»",*) +V xc {r,t) , (2) 

where the first term corresponds to the kinetic energy, 
and the following ones represent the external poten- 
tial, the Hartree potential that describes the classical 



interaction between the electrons, and the exchange- 
correlation term that accounts for all non-trivial parts of 
the electron-electron interaction. Note that the Hartree 
and exchange-correlation terms are time-dependent as 
they are functionals of the (time-dependent) density. 
This latter quantity can be evaluated from the occupied 
Kohn-Sham orbitals 



(3) 



We are concerned with external potentials that are the 
sum of a time- independent part, typically created by a set 
of nuclei, and a monochromatic electric field va c id(r, t) — 

Ei=i r i cos (ut). If we assume that the magnitude of 
A is small, we can use perturbation theory to expand 
the Kohn-Sham wave-functions in powers of A. The first 
order term reads 



ip m (r,t) = e 



{V4°V)- 



\ [Ai^SiCr, + V-^Vffi^, -w)] } , (4) 



where ipm (r) are the wave functions of the static Kohn- 
Sham Hamiltonian H^ ' obtained by taking A = 



(5) 



and tpm\( r ^ U! ) are the first order variations of the time- 
dependent Kohn-Sham wave-functions. 

From Q and the definition of the time-dependent den- 
sity Eq. ([3]), we can obtain the time dependent density 
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n (r, t) = n(°> (r) + - £ [A.e^'nf ) (r , u) 
i 

+ \ % e^ t nf\r,-u)}), (6) 
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with the definition of the first-order variation of the den- 
sity 
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4 1) (r,c)=^{[^) (r)] ^(r.w) 



^(r,-c)]*^)(r)} . (7) 



By replacing the expansion of the wave- functions (|4]) in 
the time-dependent Kohn-Sham equation fl}, and pick- 
ing up the first order terms in A, we arrive at a Stern- 
heimer equation for the variations of the wave functions 



{H^ -e m ±LJ + i V } .(r, ±to) = 

-P c H\ 1 \±w)^{r), (8) 
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with the first order variation of the Kohn-Sham Hamil- 
tonian 



h[ 1] {lo) = n+ ld 4 r 



]3 , n«(r» 



j 

- /dV/ M (r.r');, a, {r'.^). CM + i ^ - <^> w fe )AjA fe cos(w,-f) cos(o; fe i) 

+ ... (12) 
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where P c is the projector onto the unoccupied subspace. 
The effect of this projector is to make zero the compo- 
nents of ip^\ (r , ±lu) in the subspace of the occupied 
ground state wavefunctions. In linear response, these 
components do not contribute to the variation of the 
density 60 , therefore we can safely ignore the projector. 
This is important for large systems as the cost of the cal- 
culation of the projections scales quadratically with the 
number of orbit als. 

The first term of h[ 1] (w) comes from the external per- 
turbative field, while the next two represent the variation 
of the Hartree and exchange-correlation potentials. The 
exchange-correlation kernel is a functional of the ground- 
state density n^°\ and is given by the functional deriva- 
tive 



/ xc [n(°)](r,r') 



Sv xc (r) 



Sn(r') 



(10) 



In the previous equations we made use of the adiabatic 
approximation to write / xc as a frequency independent 
quantity. Equations and @ form a set of self consis- 
tent equations for linear response, that only depend on 
the occupied ground state orbitals. 

Note that we included in Eq. © a positive infinitesimal 
X). This term is essential to obtain the correct position of 
the poles of the causal response function, and therefore 
to obtain the imaginary part of the polarizability. Fur- 
thermore, using a small, but finite, r\ allows us to solve 
numerically the Sternheimer equation close to resonan- 
ces, as it removes the divergences of this equation. 

By following the same kind of reasoning, we can arrive 
at a hierarchy of Sternheimer equations for the higher 
order terms in A. These will be needed for the calculation 
of 7, the second order hyperpolarizability, or higher order 
hyperpolarizabilities.— 



B. Polarizability 

The time dependent dipolc moment is defined as 

to(t) = Jd 3 rn(r,t)n. (11) 

The polarizabilities are defined by the expansion of the 
dipole moment in terms of the electric field 



We must notice that there are several conventions for 
the definition of the (hyper Vpolarizabilities, which are 
conveniently detailed in Refl23l In this work we follow 
convention AB (where the prefactors l/n\ are explicitly 
included in Eq. [12jl . that appears to be the most used 
by the theoretical community. All referenced values have 
been converted to this convention. 

If we replace expression ([7} in Eq. (|11[) and compare 
with (|12p we can obtain a formula for the polarizability 
in terms of the variation of the density 



atij (to) = j d 3 r nf' (r , u))ri . 



(13) 



The quantity most easily accessible experimentally is 
the photoabsorption cross section, that can be evaluated 
directly from the linear polarizability 



aUS) = — 3ma(w) , (14) 
c 



where a is the trace of the polarizability tensor 



1 3 

"H = o^2 a ^) ■ ( 15 ) 
6 i=l 



C. First hyperpolarizability 

If for the dynamic hyperpolarizabilities we follow the 
same procedure as before, we get an expression in terms 
of the second order variation of the density that requires 
the evaluation of higher order variations of the wave func- 
tions. However, it is possible to get the first hyperpolari- 
zability directly from the first order variations by means 
of the 2n + 1 theorem. This theorem states that the 
nth order variations of the wave functions are enough to 
obtain the 2n + 1 derivative of the energy This the- 
orem can be expanded to the dynamic case and allows 
us to write the first hyperpolarizability (5 in terms of the 
first order variations of the wave functions. After some 
algebra, we arrive at— 
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I3 ijk (-u 1 ;^,u 3 ) = -^ E E ^r^*{T-^ x )Hf\^ 2 )^] k (r,^s) 
p C=±i I m 

occ. „ „ 



- /dV/dV/dV'i^ xc (r,r^r>f ) (r^ 1 )^^f ) (r^c. 2 )4 1) (r'^ W 3) • (16) 



r 



where the first sum is over the permutations P of the 
pairs (i,—u>i), (J, cj 2 ), and (fc,cj 3 ) and the exchange- 
correlation kernel, written in the adiabatic approxima- 
tion, reads 



K xc {r,r',r") 



=(r) 



6n(r')6n(r") 



(17) 



The first hyperpolarizability tensor has 27 components 
and is in general non-symmetric. The quantity that is 
experimentally relevant is 



1 3 



(18) 



Where z is oriented in the direction of the dipole moment 
of molecule. Sometimes the equivalent quantity /3 voc — 
f3 z = 5/3 (3\\ is used. 

III. IMPLEMENTATION 

This scheme has been implemented using a real space 
grid-based formulation in the code octopus. 25 We have 
chosen a real space grid, as it allows us a systematic con- 
vergence of the results (hyperpolarizabilities are notori- 
ously difficult to converge with localized basis sets) . How- 
ever, uniform grids can not easily describe all-electron 
atoms, so we replace the electron- nuclear Coulomb inter- 
action by Klcinman-Bylander pseudopotentials. This is, 
however, a well controlled approximation for the systems 
we are interested in. 

In (TD)DFT several approximations exist for the 
exchange-correlation termi 2 ^ In our approach we can 
treat the exchange-correlation term at two levels, one is 
the ground state exchange-correlation potential involved 
in the calculation of the ground state wave functions, 
and the other is the exchange-correlation kernel. For the 
ground state, except were noted otherwise, we use the 
local density approximation (LDA), we will also use the 
exact exchange functional in the KLI 32 approximation. 
For the exchange-correlation kernel, we have decided to 
use, due to its simplicity, the adiabatic local density ap- 
proximation (LDA). (Although our scheme is quite gen- 
eral and can in principle be applied to any exchange- 
correlation functional.) 



The LDA is a well studied approximation, and is quite 
reliable in the prediction of many properties. One impor- 
tant main defect in this context is the wrong asymptotic 
part of the LDA exchange-correlation potential, that for 
neutral systems decays exponentially instead of falling 
as 1/r. This usually leads to small HOMO-LUMO gaps 
which implies systems that are too polarizable, in con- 
trast to Hartree-Fock where the gap is larger and the 
magnitudes of the polarizabilities are underestimated. 
The exchange and correlation kernel will contribute to 
reduce the independent particle polarizability even at 
ALDA level (this contribution for extended systems is 
zero and as the long range behavior of the XC poten- 
tial is not relevant in this regime, this clearly points 
that the non-locality of the XC kernel as well as self- 
interaction correction are responsible for the bad perfor- 
mance of LDA). We will observe this overestimation in 
the calculations that follow. In fact, in the case of com- 
pact finite systems there has been indications that the 
exchange-correlation potential seems to be more impor- 
tant than the kernel^ The situation is particularly prob- 
lematic in the case of long molecular chains^ 2 - 7 - where stan- 
dard exchange-correlation functionals can greatly over- 
estimate polarizabilities when compared to many-body 
approaches. Note, however, that this is not a deficiency 
of DFT, but of the LDA approximation (and of many 
exchange-correlation functionals), that can in principle 
be treated 2 ^ by using more sophisticated orbital depen- 
dent functionals like the self-interaction corrected LDA 30 
or the exact exchanged However in present orbital func- 
tionals there is still a significant correlation contribution 
that is not taken into account and is responsible for the 
discrepancies between theory and experiment for long- 
chains in the exact exchange approach. 27 

Numerically, the central part of our scheme is the solu- 
tion of the Sternheimer equation ([5]) . This has the form of 
a linear equation were the operator to invert is the shifted 
ground state Hamiltonian. As the shift is complex (due 
to the ir\ term), this operator is not Hermitian. There- 
fore, we can not use standard techniques common in the 
community, like the simple conjugated gradients scheme, 
but have to rely on more general (and involved) linear 
solvers. Our choice was the biconjugate gradient stabi- 
lized method^ 3 - Close to the resonance frequencies, the 
Sternheimer equation becomes very badly conditioned, 
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and the solution process turns out to be very costly. The 
problem can be eased by the use of preconditioning. We 
have found that a smoothing preconditioner— can dra- 
matically improve convergence in these cases. Also for 
small system the solution process can be made less costly 
if we solve the Sternheimer equation in the space of the 
unoccupied wave functions by orthogonalizing the right- 
hand side of the Sternheimer equation with respect to 
the occupied wave functions. Although this would not be 
practical for large systems as the orthogonalization pro- 
cess can become very demanding. Even with these tech- 
niques, the process is much more costly for frequencies 
near resonance. For example, for the CO case, the full 
self-consistent solution of the Sternheimer equation for a 
single frequency in a single direction requires around 1700 
applications of the Hamiltonian, for the a near resonance 
frequency approximately 10 times more applications are 
required. As a comparison, for the ground state DFT cal- 
culation around 2000 Hamiltonian operations are needed. 

As the right-hand side of the Sternheimer equation de- 
pends on the linear variation of the density, the problem 
has to be solved self-consistently. For this, we use similar 
strategies as for the ground-state calculation, mixing the 
linear variation of the density using a Broyden scheme^ 
in order to speed up the convergence of the self-consistent 
cycle. 

The Poisson equation is solved using the interpolating 
scaling functions scheme proposed in Refl36l As this pro- 
cess has to be done only once per self-consistency itera- 
tion, the performance of the Poisson solver is not critical. 

The total cost of calculating the response is of or- 
der 0(Nk s NgM UJ ), where Nk s is the number of Kohn- 
Sham orbitals, N g is the number of grid points and M u 
is the number of frequencies we desire (which is inde- 
pendent of the system size). This scaling is much bet- 
ter than for the approaches that relay on expansions in 
particle-hole states, and even better than for the ground- 
state calculation that normally scales as 0(Nk s Ng) (due 
to the necessity of orthogonalizing the wave- functions) . 
We believe, therefore, that this method can be used 
to study (hyper)polarizabilities of very large systems, 
like nanostructurcs or molecules with biological inter- 
est. After obtaining the linear response, the evaluation of 
the hyperpolarizability from Eq. (|16p has a cost propor- 
tional to 0(N^ s N g M UJ ) but with a very small prefactor. 
Note that O(N) (with N the number of atoms) schemes 
are available for ground stated! and static polarizability 
calculations^ These linear scaling methods are based 
on "nearsightness"— or the idea of divide and conquer 4^ 
The real space treatment used in our method would allow 
for the incorporation of these strategies to reduce even 
further the numerical cost for large systems. 



IV. EXAMPLES OF APPLICATIONS 

We decided to illustrate the implementation of our 
method by applying it to some simple molecules that 
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TABLE I: Comparison of static polarizabilities and hyper- 
polarizabilities for CO. Results are in atomic units. "Finite 
differences results from Refl44l. "Experimental result from 
Refl45l. C LDA basis set results from Refill. "Experimental 
result from Refl47l. e LDA basis set results from Refl48l. 

have been well studied both experimentally and theore- 
tically: CO, H2O, and paranitroaniline (PNA). 

In all our calculations we took the experimental geo- 
metries: CO = 1.13 A; HO = 0.957 A, HOH = 104.5°; 
the experimental geometry of PNA determined through 
X-ray crystallography can be found in Ref[4l|. However, 
there are two CH distances missing from the crystallo- 
graphic data. For these we followed Refill and took their 
theoretical values (see Table 1 of Refl43l) . In all cases the 
dipole moment of the molecule was taken perpendicular 
to the z axis, for H2O the molecule was considered in the 
yz plane and in the case of PNA it was taken in the xz 
plane. We used Troullier-Martins norm-conserving pseu- 
dopotentials with a core radius of 0.66 A for H, 0.78 A 
for C and 0.74 A for both N and O. 

We used a simulation box composed of spheres cen- 
tered at each atomic position, with a radius of 9.5 A for 
CO, 7.4 A H 2 and 5.3 A for PNA. The points were dis- 
tributed in a regular grid with spacing 0.20 A for CO, 
0.17 A for H 2 0, and 0.19 A for PNA. With these para- 
meters, hyperpolarizabilities are converged to better than 
1%. As expected, the simulation boxes required to con- 
verge the hyperpolarizabilities were much larger than the 
ones typically used in ground-state calculations, as these 
quantities have sizeable contributions from the regions 
far away from the nuclei. The required simulation box 
size depends on the frequency of the perturbation, be- 
ing larger for higher frequencies. Finally, we used the 
LDA parameterization by Perdew and Zunger— to ap- 
proximate the exchange-correlation functionals. 

We will start our discussion with CO. The results for 
static properties are displayed in Table HI The results 
fully agree with the other DFT-LDA calculations. Com- 
pared to more sophisticated quantum chemistry meth- 
ods, the LDA overestimates the values for the polarizabil- 
ities and the hyperpolarizabilities. This, as mentioned 
before, comes from a deficiency in the asymptotic region 
of the LDA potential. 

We now turn to the dynamic properties. In Fig. [1] 
we plot the absorption cross section of the CO molecule 
obtained with our approach, and compared to the spec- 
trum obtained through the direct solution of the time- 
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FIG. 1: Average photo-absorption cross section of the CO 
molecule, calculated within the adiabatic LDA. The line cor- 
responds to the results obtained through the solution of the 
Sternheimer equation, while the dots are obtained through the 
solution in real-time of the time-dependent Kohn-Sham equa- 
tion. Low resolution experimental absorption cross-section 
from Refl49l. For a more detailed comparison, the relevant ex- 
perimental excitation energies are at: 8.51 eV (A 1 !!), 10.78 eV 
(B 1 ^), 11.40 eV (C 1 S+), and 11.53 eV (E 1 ^^- 



dependent Kohn-Sham equations in real-time. As ex- 
pected, the two calculations agree perfectly, which vali- 
dates our numerical implementation of the Sternheimer 
equations. Note that, even if both methods have identi- 
cal scalings, the solution of the Sternheimer equation still 
has a larger prefactor if the whole spectrum is required. 
This is due to the ill-conditioning of the linear system 
close to resonances. From Fig. [1] it is clear that the two 
theoretical results are red-shifted with respect to the ex- 
perimental curve. This shortcoming of the simple LDA 
can be corrected by using functionals with the correct 
asymptotic behavior. 29 

In Fig. [2] we present our calculations of the second har- 
monic generation spectrum, /3u(— 2lo;lo,lo), of CO, to- 
gether with the available experimental results^ and pre- 
vious theoretical data42- Our results agree very well with 
previous DFT results using the LDA. We can also see that 
the use of the generalized gradient approximation BLYP 
(Becke 88 51 for exchange and Lee, Yang, and Parr- 
correlation) does not change significantly the results. - 
Using a hybrid functional, the Becke 3 parameter B3LYP 
functional^ does reduce the error, while Hartree-Fock 
(HF) results underestimate the value for the hyperpo- 
larizability. The best results are, as expected, obtained 
by coupled cluster calculations using singles and doubles 
(CCSD). 

Now we turn to the H2O molecule. To test our im- 
plementation, we show, in Table HH the different com- 
ponents of the hyperpolarizability tensor for lu = OeV, 
1.79 eV, and 1.96 eV. We see that our results compare 
well to previous theoretical work42 The small difference 



FIG. 2: (Color online) Second harmonic generation 
/3i|(— 2uj; tjj,u)) of CO. The inset shows a comparison of the 
results of this work (TDLDA) with other available results. 
Exp: Experimental results from Ref@; HF: Hartree Fock re- 
sults from Ref|H; LDA, BLYP, B3LYP: DFT results from 
Ref|H CCSD: Coupled cluster results from RefS2. 
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TABLE II: Tensor components for second harmonic genera- 
tion /3(-2w;w,w) of H 2 0. "Theoretical results from Ref|H. 



can be explained by the different numerical methodolo- 
gies (real-space grid and pseudopotentials in our case and 
basis sets in Refl42l). 

In Table |TTI] we show the static polarizability, a(0, 0), 
the static first hyperpolarizability, j3» (0; 0, 0), the optical 
rectification, /3ii(0;cj, —lu), and second harmonic genera- 
tion, /3n (— 2lu; lu, lu) for water. All dynamic values were 
calculated for lu = 1.79 eV. This time we also have in- 
cluded results with the KLI orbital dependent exchange- 
correlation potential combined with the ALDA exchange- 
correlation kernel. Concerning LDA values, we can see 
that our results fully agree with the results of Reffl7l that 
also uses a grid based representation, while basis-set re- 
sults from Ref.42 differ by less than 10%. We can see 
that all LDA and BLYP values overestimate the mag- 
nitude of (hypcr)polarizabilities with respect to experi- 
ment and coupled cluster calculations. The use of more 
sophisticated exchange correlation functionals, like LB94 
or B3LYP, improves significantly the results, while the 
KLI/ ALDA scheme gives an underestimation of the mag- 
nitude of (hypcr)polarizabilities similar to Hartee-Fock 
results. 
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q(0,0) /3||(0;0,0) ^(Ojw.-w) /? n (- 


-Zuj, uj, oj) 


This work (TDLDA) 


10.51 


-25.89 


-28.33 


-34.71 


This work (KLI/ALDA) 


8.61 


-11.75 


-12.43 


-14.03 


LDA™ 


10.5 


-26.1 


-28.6 


-35.1 


BLYP a 


10.8 


-27.9 


-30.9 


-38.8 


LB94 a 


9.64 


-17.8 


-17.7 


-20.3 


LDA b 


10.63 


-23.78 


-26.09 


-32.12 


BLYP 6 


10.77 


-23.65 


-26.11 


-32.76 


B3LYP b 


9.81 


-18.54 


-20.11 


-24.11 


HF b 


8.53 


-10.73 


-11.27 


-12.52 


CCSD(T) C 


9.79 


-18.0 


-19.0 


-21.1 


exp. 


9.81 d 






-22±9 e 



TABLE III: Comparison of (hyper)polarizabilities of H2O for different DFT calculations. CCSD(T) and experimental values 
are Riven as reference. For dynamic results w = 1.79 [eV]. "Grid based calculations from Refil7l; 6 Basis set calculations from 
Refl42l: c Results from Refl53: Experimental results from ReflHEl: Experimental results from Refl56l: 
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FIG. 3: Calculated optical rectification /3\\ (0;a>; —oj) of H2O. 



To conclude the discussion of water results, we plot, in 
Fig. [3] the frequency dependence of the optical rectifica- 
tion in the visible and near ultraviolet regimes. It is clear 
that our method not only works for small non-resonant 
frequencies but also in the more complicated resonant 
regime. 

Finally, we turn to a larger molecule: paranitroaniline. 
Our results for the second harmonic generation process 
in this molecule are given in Fig. [4j There are two expe- 
riments available: i) a gas phase experiment^ performed 
for a single frequency; and ii) paranitroaniline in solvent 
for several frequencies M- The latter were corrected for 
the presence of the solvent, but this correction is clearly 
incomplete as the value from Ref 57 for oj = 1.17eV is 
still 10% larger that the gas phase measurement^ We 
include also other theoretical values using DF T 48 i 59 and 
CCSD.— We can observe that our results underestimate 
the solvent experimental results by about 15% for all 
available frequencies. In comparison, B3LYP and CCSD 



results are seriously too small at high frequencies, with 
values that are around 40-50% smaller than experiment. 
We think that the reason for this discrepancy is the bet- 
ter description of the hyperpolarizabilities near resonance 
of our method. Furthermore, it uses a grid based repre- 
sentation that describes better the regions far from the 
nuclei in comparison with localized basis sets, allowing 
for larger flexibility in capturing the dynamic changes in 
the wave functions. 



V. CONCLUSIONS 

In summary, we have presented a method that allows 
the calculation of both static and dynamic polarizabili- 
ties and hyperpolarizabilities. Our approach is based on 
the Sterheimer equation, within the formalism of time- 
dependent density functional theory, and requires the so- 
lution of an non-Hermitian linear equation. This solution 
is obtained through a generalization of the conjugated 
gradients method using a real-space (basis set free) rep- 
resentation of the wave- functions. In this way we are able 
not only to obtain static quantities, but also the whole 
frequency dependence of the (hyper)polarizabilities even 
close to resonances. The scaling with the number of 
atoms in the system is excellent, so we expect that the 
method will be useful for the study of very large sys- 
tems. First applications to small benchmark molecules 
yield quite good results in comparison to previous theo- 
retical approaches and experimental results. 

One of the beauties of this approach is how easily it can 
be generalized to higher orders and to handle other kinds 
of static or dynamic perturbations. For example, phonon 
frequencies, (resonant) Raman tensors, NMR tensors, 
forces in the excited state, etc. can all be obtained by just 
changing the right-hand side of the Sternheimer equa- 
tions. The third and higher order polarizabilities can 
also be obtained by solving a hierarchy of Sternheimer 
equations that have the same form as the first order one. 
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FIG. 4: (Color online) Second harmonic generation f3\\ (— 2uj; uj, uS) of paranitroaniline. Note that the y axis in the right panel 
is in logarithmic scale, exp. solv.: Solvent phase experimental results from Refl57l; exp. gas: Gas phase experimental results 
from ReflU; LDA/ALDA, LB94/ALDA: DFT basis set results from Ref|H; B3LYP: DFT basis set results from Ref|H; CCSD: 
Coupled cluster results from Refl42l. Some references use a different convention to define hyperpolarizabilities, all the values 
shown here have been converted to convention AB. 



Work has already started to extend our implementation 
in these directions. 
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